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Source  localization  in  shallow  water  using 
multi-frequency  processing  of  shot  data 

G.  Haralabus  and  P.  Gerstoft 


Executive  Summary:  Single-frequency  techniques  which  proved  to  be  successful 
for  source  localization  in  deep-water  are  shown  to  be  inadequate  in  littoral  areas.  To 
overcome  this  problem,  the  processing  schemes  need  to  be  redesigned  to  extract  and 
incorporate  more  information  from  the  received  signal  and  improve  the  estimates  of 
parameters  which  describe  the  propagation  channel.  This  report  describes  the  use  of 
multi-frequency  processing  methods  in  a  global  search  scheme  in  which  the  source 
coordinates  are  estimated  jointly  with  environmental  parameters.  The  analysis  is 
based  on  comparisons  between  modelled  data  and  real  data  from  explosive  sources 
collected  during  a  SACLANTCEN  sea  trial  in  the  Mediterranean  Sea,  in  October 
1993. 

The  variability  of  geometric  (source  coordinates,  array  tilt),  and  environmental  (water 
depth  at  the  source  and  the  receiver),  parameters  as  a  function  of  the  number  of 
frequencies  is  examined.  It  is  shown  that,  for  the  particular  data  set,  stable  and 
reliable  results  are  obtained  when  at  least  ten  frequencies  with  the  highest  energy 
content  are  used.  Beyond  this  ‘saturation’  point  the  estimates  remain  relatively 
constant.  Based  on  these  observations,  the  localization  performance  of  a  multi¬ 
frequency  version  of  the  Bartlett  processor  is  examined.  It  is  demonstrated  that 
insufficient  knowledge  of  the  propagation  channel  limits  the  resolution  of  the  source 
location  for  both  single-frequency  and  multi-frequency  cases.  The  use  of  global 
inversion  methods  based  on  genetic  algorithms  for  the  joint  optimization  of 
environmental  and  geometric  parameters  is  shown  to  considerably  improve  the 
localization  performance  of  the  processor. 

The  suggested  multi-frequency  methodology  is  related  primarily  to  the  way  the 
received  signal  is  processed  and  is  independent  of  the  geometry  of  the  problem, 
therefore  it  can  be  equally  well  applied  to  both  active  and  passive  sonar. 
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Source  localization  in  shallow  water  using 
multi-frequency  processing  of  shot  data 

G.  Haralabus  and  P.  Gerstoft 


Abstract:  Multi-frequency  processing  methods  are  applied  to  real  data  generated 
by  explosives  to  examine  the  variability  of  geometric  and  environmental  parameter 
estimates.  The  frequencies  are  combined  in  an  incoherent  fashion  and  genetic 
algorithms  are  employed  in  an  efficient  global  search  scheme. 

The  source  coordinates,  the  water  depth  at  the  source  and  the  receiver  location,  and  the 
tilt  of  the  vertical  array  are  estimated  as  functions  of  the  number  of  frequencies 
processed.  Using  three  different  processing  schemes,  it  is  found  that  when  less  than 
ten  frequencies  are  employed  the  estimates  are  unstable,  whereas  after  this  critical 
number  a  constant  estimation  level  is  attained. 

The  localization  performance  of  a  multi-frequency  version  of  the  Bartlett  processor  is 
also  examined.  The  degradation  of  the  results  due  to  a  non-optimized  environment  is 
demonstrated.  The  source  resolution  is  shown  to  be  improved  when  the  search  for  the 
source  coordinates  incorporates  simultaneous  estimation  of  the  environmental 
parameters. 
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1 

Introduction 


Acoustic  field  inversion  using  matched  field  processing  has  been  successfully  used 
to  determine  both  geoacoustic  and  environmental  parameters.  It  has  been  found 
that  the  localization  performance  increases  when  selected  waveguide  parameters  are 
determined  simultaneously  [1],  It  has  been  shown  that  the  stability  of  the  parameter 
estimates  is  significantly  improved  when  using  data  at  a  few  well  defined  frequencies 
[2],  [3].  A  disadvantage  of  this  approach  is  that  the  computational  time  increases 
linearly  with  the  number  of  frequencies  used  in  the  inversion  procedure.  This  study 
examines  multi-frequency  processing  strategies  for  a  transient  signal  in  a  mildly 
range-dependent,  shallow  water  environment. 

The  data  for  this  study  were  collected  during  a  SACLANTCEN  sea  trial  in  1993  [1]. 
More  information  about  the  experiment  is  given  in  the  next  section.  Replica  fields 
are  calculated  using  the  SNAP  [4]  adiabatic  normal  mode  program.  The  environment 
can  be  described  as  a  half-space  divided  into  a  water  column,  a  sediment  layer,  and 
a  semi-infinite  sub-bottom  (Fig.  1).  It  can  be  seen  that  this  is  a  range  dependent 
scenario  in  which  the  bathymetric  difference  between  the  source  and  the  receiver 
is  approximately  7  m.  Comparisons  between  real  and  synthetic  data  are  based  on 
the  Matched  Field  Processing  (MFP)  technique  [5],  [6].  This  method  compares 
the  relative  phase  and  amphtude  of  recorded  signals  with  model  predictions.  To 
find  the  best  fit  between  real  and  synthetic  data,  modelling  and  matching  tools  are 
incorporated  into  the  SAGA  [7]  global  inversion  algorithm;  a  time-efficient,  multi¬ 
dimensional  search  method  based  on  genetic  algorithms.  In  general,  the  following 
procedure  is  followed:  a)  The  bounds  of  the  parameters  to  be  estimated  are  set  based 
on  a  priori  knowledge.  Then,  the  environment  is  discretized  and  fed  into  SNAP,  b) 
SNAP  calculates  the  acoustic  field  at  the  receiver  location,  c)  Real  and  simulation 
data  (transformed  in  the  frequency  domain)  are  compared  by  means  of  an  objective 
function,  d)  The  above  steps  are  repeated  for  the  next  chosen  parameter  set.  The 
iteration  procedure  is  controlled  by  SAGA,  e)  Statistical  analysis  and  a  posteriori 
probability  plots  of  the  obtained  estimates  are  provided. 
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Source  Receiver 


Figure  1:  Measured  sound-speed  profile,  bathymetry,  and  geoacoustic  parameters 
for  the  experimental  site  north  of  the  island  of  Elba. 
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Data  collection 


The  data  were  collected  during  a  SACLANTCEN  sea  trial  north  of  the  island  of  Elba 
in  the  Mediterranean  Sea,  in  October  1993.  Twenty-five  explosive  sound  sources 
(SUS  MK-64,  model  0)  [8],  [9]  were  deployed  during  a  triangular  course  with  sides 
of  23  km.  The  propagated  signals  were  received  by  a  vertical  array  which  was  sit¬ 
uated  at  the  eastern  apex  of  the  triangle.  Shot  5  was  singled  out  for  examination 
because  this  particular  explosion  rendered  a  better  recording  than  the  others,  and  is 
also  investigated  by  Tolstoy  et  al  [10],  using  different  methods  than  those  described 
in  this  paper.  The  number  of  hydrophones  used  and  the  sampling  rate,  were  deter¬ 
mined  by  the  equipment  on  board  the  research  vessel.  For  such  practical  reasons, 
in  this  experiment,  the  sampling  frequency  was  1000  Hz.  Therefore,  the  maximum 
frequency  is  500  Hz  although  it  is  known  that  the  frequency  band  of  explosions  ex¬ 
tends  up  to  10  kHz  [11],  [12].  Figure  2  shows  the  time  series  and  the  spectrum  of  the 
signal  received  on  32  hydrophones  (with  2  m  spacing)  across  the  vertical  array.  The 
explosives  were  deployed  from  a  ship  which  followed  a  preassigned  path,  so  there  is 
probably  an  offset  between  the  registered  coordinates  and  the  site  of  the  explosion. 
It  would  have  been  useful  to  have  had  additional  sound  velocity  profile  and  bathy¬ 
metric  measurements  along  the  propagation  path  because  acoustic  propagation  is 
strongly  dependent  on  these  factors.  Table  1  shows  the  five  parameters  considered  in 
this  study,  the  search  interval  for  each,  and  the  baseline  values  for  this  experiment. 


Parameter 

S.  range 

5.  depth 

Depth  at  S. 

Depth  at  R. 

Tilt 

Interval 

11-13  km 

10-20  m 

120-128  m 

114-121  m 

-3-3  m 

Baseline 

12.1  km 

18  m 

122  m 

115  m 

0  m 

Table  1:  Environmental,  source  and  receiver  parameters,  and  baseline  model  values. 
S  and  R  denote  source  and  receiver  respectively. 
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Figure  2:  Time  series  (a)  and  frequency  (b)  spectrum  of  the  signal  received  on  the 
vertical  array. 
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Objective  function 


Parameter  estimation  is  treated  as  an  optimization  problem  in  which  the  algorithm 
searches  for  the  model  vector  m  which  minimizes  an  objective  function  $(m).  For 
a  single  frequency  u,  let  d(a;)  be  the  real  data  received  on  a  vertical  array  of  N 
hydrophones,  and  p(u;)  the  data  predicted  from  the  simulation  model.  Real  and 
synthetic  data  are  related  by  the  equation 


d(w)  =  p(w)  +  n(w)  (1) 

where  n(a;)  represents  the  noise.  Since  the  synthetic  data  can  be  expressed  as 

p(w)  =  w(<.u,  m)5(a;)  (2) 

where  w(a;,  m)  is  the  transfer  function  of  the  acoustic  medium  and  5(w)  is  the 
deterministic  source  spectrum,  Eq.  (1)  becomes 


d(ci;)  =  w(u;,  m)5(u;)  +  n(u;) 


(3) 


or 


n(a;)  =  d(a>)  —  w(cu,  m)s(u;) 


(4) 


The  n  term  includes  primarily  environmental  noise  and  other  factors  of  uncertainty, 
such  as  modelling  and  measurement  inaccuracies.  These  factors  of  error  are  consid¬ 
ered  to  be  independent,  random  processes  without  burst-like  quality.  It  is  therefore 
reasonable  to  assume  that  n  has  a  Gaussian  distribution  with  zero  mean  and  covari¬ 
ance  matrix  C,  i.e.,  n  G  N(0,C).  Consequently,  the  real  data  d  are  also  Gaussian 
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distributed,  with  mean  ws  and  the  same  covariance  matrix  C,  i.e.,  d  G  N(w5,C) 
(for  simplicity  the  dependence  on  lj  and  m  is  omitted). 

For  complex  data  the  probability  density  function  is  defined  as: 


p(d)  = 


— (d  —  ws)^C  ^(d  — ws) 


(5) 


where  N  is  the  vector  dimension,  f  is  the  complex-conjugate  operator,  and  |C| 
denotes  the  deternimant  of  C  matrix. 

Assuming  that  the  noise  is  spatially  uncorrelated  and  an  identical  spectrum  for  each 
hydrophone,  the  cross-spectral  covariance  matrix  can  be  written  as  C  =  ul,  where 
the  noise  standard  deviation  n  depends  only  on  frequency  and  I  is  the  identity 
matrix. 

In  the  case  of  broadband  data,  the  processor  incorporates  individual  frequencies  in  an 
incoherent  fashion.  Given  the  observed  data  sets  d;  for  different  frequencies  w/  and 
assuming  that  they  are  uncorrelated  with  each  other,  the  optimization  problem  is 
modified  in  order  to  find  the  model  vector  m  which  creates  a  replica  field  w  =  w(m) 
which  minimizes  the  logarithm  of  the  likelihood  function 


M  =  ^/o5(p(d;))  = 

1=1  1=1 


—  Nlog-KUi  -  — (d;  —  W;S/)^(d;  —  W;S;) 
Vl 


(6) 


where  I  =  1,...,X  denotes  the  number  of  frequencies.  When  the  expression  for  a 
harmonic  point  source  si  —  is  substituted  in  Eq.  (6)  and  one  solves  the 

equations  =  0  and  =  0,  the  estimate  of  the  source  signal  becomes 


w/d/ 

W^W/ 


(7) 


By  substitution  of  Eq.  (7)  into  Eq.  (6),  the  likelihood  function  is  modified  as: 
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M=i: 


-Nlogiri/i  — 


VI 


d/d, 


w|d,d/w, 

wjwi 


(8) 


For  each  frequency  the  cross-spectral  covariance  matrix  of  the  actual  data  is  calcu¬ 
lated  from  the  formula 


R,  =  d/dj 


(9) 


which  implies  that  its  trace  is 


trR,  =  d/d, 


(10) 


Using  Eq.  (9)  and  (10),  Eq.  (8)  becomes 


L 


m  =  y, 

l=i 


-Nlog-Kui - 

VI 


w/R,w, 

w/w. 


(11) 


This  expression  can  be  further  simplified  depending  on  the  type  of  assumption  which 
is  made  for  the  noise  v: 

A)  Assuming  that  u  is  not  constant,  then  by  solving  the  equation  =  0,  it  is 
derived  that 


w1r,w, 

w/w. 


(12) 


When  Eq.  (12)  is  substituted  in  Eq.  (11)  and  unnecessary  constants  are  removed, 
the  exponential  form  of  the  likelihood  operator  is  defined  as  the  problem’s  objective 
function 
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L 


$=n 

/=i 


trR/  — 


w|R/W/ 

w/w; 


(13) 


B)  Alternatively,  if  v  is  assumed  constant,  then  directly  from  Eq.  (11),  the  objective 
function  becomes 


L 


;=i 


trR;  — 


wJR;W; 

wjw; 


(14) 


In  both  cases,  the  algorithm  is  searching  for  the  model  vector  m  which  generates 
a  synthetic  field  w(m)  minimizing  #(m).  In  the  present  work  we  utilize  Eq.  (14). 
Notice  that  the  second  term  in  Eq.  (13)  and  (14)  is  the  expression  for  the  Bartlett 
processor  [5],  i.e.,  for  each  frequency  ui 


B;  =  wJR/W;  = 


N 

E 

i=l 


(15) 


It  should  be  noted  that  the  same  formulae  apply  when  the  real  data  used  are  created 
from  the  average  of  several  time  segments.  Then,  the  covariance  matrix  can  be 
calculated  according  to  the  formula 


R,  =  4^d,d*  (16) 

where  K  is  the  number  of  time  segments. 

Finally  it  should  be  mentioned  that  for  each  frequency,  the  objective  function  is 
based  on  the  relative  phase  differences  (weighted  with  pressure  magnitude)  across 
the  elements  of  the  receiving  array.  The  L  frequencies  used  in  the  inversion  should 
be  the  ones  that  are  best  received  on  most  hydrophones.  To  satisfy  this  condition, 
we  calculate  the  array  averaged  spectrum  in  which  each  frequency  bin  corresponds 
to  the  summation  of  the  energy  components  from  aU  32  hydrophones.  Based  on  this 
calculation,  the  frequencies  corresponding  to  the  maximum  values  are  selected. 
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Single-frequency  parameter  estimation 


To  demonstrate  the  need  for  multi-frequency  analysis,  one  must  provide  evidence 
that  single  frequency  algorithms  are  insufficient  for  accurate  parameter  estimates. 
As  shown  in  this  section,  methods  based  on  a  single  frequency  lead  to  unstable 
results. 

Here,  two  single-frequency  scenarios  for  146.5  and  194.3  Hz  are  examined.  First  the 
Fourier  transform  of  the  complete  signal  is  calculated  and  then  for  each  hydrophone 
the  values  corresponding  to  146.5  and  194.3  Hz  are  selected.  For  each  frequency,  the 
covariance  matrices  are  calculated  according  to  Eq.  (9).  In  both  cases  the  search 
bounds  and  the  baseline  model  are  the  same.  The  SAGA  output  is  shown  in  Fig. 
3.  The  vertical  lines  represent  the  baseline  model  and  the  a  posteriori  probability 
plots  indicate  the  best  parameter  estimates. 

It  is  obvious  that  the  two  predictions  are  quite  different.  Also,  it  is  important  to 
notice  that  the  a  posteriori  probabilities  have  low  standard  deviation.  This  signifies 
that  the  algorithm  did  converge  with  the  best  solution  for  each  frequency;  in  other 
words  an  increase  in  the  number  of  iterations  would  not  change  the  accuracy  of 
these  predictions.  Key  reasons  for  having  inconsistent  parameter  predictions  when 
different  frequencies  are  processed  are  the  insufficient  knowledge  of  the  propagation 
waveguide  (e.g.  sparsely  collected  environmental  measurements),  and  the  small 
number  of  time  samples.  The  next  step  is  to  obtain  more  consistent  estimates  by 
combining  information  across  the  entire  frequency  spectrum. 
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Figure  3:  A  posteriori  distribution  based  on  one  frequency:  a)  146.5  Hz,  b)  194.3 
Hz. 
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Multi-frequency  parameter  estimation 


To  overcome  the  inaccuracy  of  single-frequency  estimates  shown  in  the  previous  sec¬ 
tion,  multi-frequency  scenarios  are  investigated.  The  global  inversion  algorithm  is 
now  based  on  different  data  sets  which  include  one  to  twenty  five  frequencies.  These 
frequencies  are  selected  according  to  energy  content.  Multi-frequency  processing 
is  performed  in  an  incoherent  fashion  as  indicated  by  Eq.  (14).  The  object  is  to 
obtain  estimates  which  are  close  to  the  baseline  model  and  remain  unchanged  be¬ 
yond  a  “saturation”  point.  These  two  conditions  ensure  reliable  results  and  reduces 
computation  time. 

Initially,  the  calculation  of  cross-spectral  covariance  matrix  is  based  on  the  Fourier 
transform  of  the  complete  set  of  time  samples.  For  each  frequency,  the  spectrum 
values  are  calculated  for  every  hydrophone  and  from  Eq.  (9)  the  covariance  matrix 
is  estimated.  This  procedure  is  repeated  for  all  frequencies  and  then  the  parameters 
are  estimated  using  Eq.  (14).  The  results  from  this  process  are  shown  in  Fig.  4.  The 
source  coordinates  correspond  to  the  first  two  plots.  Notice  that  when  the  number 
of  frequencies  is  less  than  10,  the  predictions  for  the  source  range  and  depth  are 
unstable.  Beyond  this  point  the  estimates  reach  a  plateau  which  is  close  to  baseline 
value.  The  predictions  for  the  water  depth  at  the  source  and  the  receiver  do  not 
follow  the  same  trend.  The  first  is  overestimated  for  aU  frequency  numbers  while 
the  second  has  sudden  fluctuations  around  the  critical  point  of  10  frequencies.  Later 
it  is  shown  that  this  fluctuation  problem  is  eliminated  by  using  different  processing 
schemes.  Finally,  the  array  tilt  demonstrates  smooth  convergence  toward  its  mean 
position. 

Obviously,  the  multi-frequency  results  are  superior  to  the  single-frequency  ones  be¬ 
cause  they  are  in  general,  stable  and  accurate.  As  expected,  when  the  number  of 
frequencies  increases,  the  predictions  gradually  improve  until  they  reach  a  stable 
value.  This  behaviour  is  intuitively  more  robust  whether  the  estimates  are  close  to 
the  baseline  model  (cases  of  source  coordinates  and  array  tilt),  or  are  not  close  to  it 
(case  of  water  depth  at  the  source  depth). 

This  analysis  is  based  on  the  spectrum  shown  in  Fig.  2.  The  stability  of  the  estimates 
depends  on  the  time  variability  of  this  energy  distribution.  To  demonstrate  this 
variability,  the  array  averaged  spectrogram  [13]  is  calculated  for  a  256  ms  time 
window  which  slides  over  the  time  domain  with  a  10  ms  increment  (Fig  5).  It 
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Figure  4:  Parameter  estimates  as  a  function  of  the  number  of  frequencies.  Unmarked 
solid  lines  indicate  the  baseline  model  values.  The  whole  time  series  is  used  and  no 
averaging  is  performed. 
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Figure  5:  Array-averaged  spectrogram. 
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can  be  seen  that  the  energy  content  remains  at  the  same  level  for  only  a  short 
period  (from  250  to  450  ms)  corresponding  to  the  high  energy  segment  of  the  signal. 
Parameter  estimates  based  on  this  time  segment  are  expected  to  be  very  stable.  To 
take  advantage  of  this  idea,  two  alternatives  can  be  proposed:  a)  to  process  only  the 
high  energy  segment  of  the  returned  signal  between  the  250th  and  the  450th  time 
sample,  and  b)  to  use  an  averaging  technique  based  on  the  entire  time  series.  In 
both  cases,  the  main  idea  is  to  eliminate  or  reduce  the  effect  of  the  “noise-like”  low- 
energy  segment  of  the  signal  between  the  450th  and  the  650th  sample,  and  increase 
the  influence  of  the  strong  signal. 

In  the  first  alternative,  the  estimates  of  the  five  parameters  of  interest  are  presented 
in  Fig.  6.  The  source  coordinates  are  shown  in  the  first  two  plots.  Similarly  to  the 
previous  scenario,  for  less  than  10  frequencies  the  estimates  are  unstable  but  beyond 
this  critical  point  both  parameters  reach  a  plateau  with  small  variations  around 
the  actual  values.  The  most  problematic  parameter  remains  the  water  depth  at 
the  source  location  which  is  overestimated.  The  water  depth  at  the  receiver  has 
a  highly  accurate  estimate  and  the  fluctuations  observed  in  the  previous  method 
is  eliminated.  Finally,  the  array  tilt  slowly  converges  with  its  mean  value  zero. 
Comparing  the  two  multi-frequency  approaches,  it  is  concluded  that  for  this  data 
set,  exploiting  only  the  high  energy  part  of  the  shot  is  the  preferred  method  as  the 
results  are  more  stable  and  the  behavior  of  the  processor  is  more  predictable. 

Proceeding  with  the  second  alternative,  the  information  from  both  segments  is  com¬ 
bined  using  averaging  techniques  to  balance  the  contribution  of  the  low  energy  part 
of  the  signal.  Here,  the  samples  between  201  and  700  are  selected.  This  series  is 
divided  into  five  non-overlapping  segments  using  a  100-sample  rectangular  window. 
Covariance  matrices  are  calculated  for  each  segment  and  each  frequency.  The  cor¬ 
responding  matrices  are  averaged  to  derive  the  final  form  of  R(w;).  The  parameter 
estimation  results  are  shown  in  Fig.  7.  The  critical  number  of  frequencies  remains  at 
ten.  It  can  be  seen  that  the  estimates  for  the  source  coordinates  are  very  close  to  the 
expected  values.  The  water  depth  at  the  source  location  is  overestimated,  as  with 
the  two  previous  cases.  The  estimate  of  the  water  depth  at  the  receiver  is  exact  and 
remains  constant  in  the  plateau  zone.  Finally,  the  array  tilt,  although  it  remains 
close  to  the  nominal  value,  does  not  have  as  good  an  estimate  as  in  the  first  method. 
In  conclusion,  the  main  characteristics  of  this  approach  are:  a)  transition  regions 
which  remain  the  same  for  aU  parameters,  and  b)  plateau  areas  which  are  flatter 
than  the  ones  in  the  previous  cases.  Both  alternatives  have  comparable  performance 
and  are  preferable  to  the  original  approach. 
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Figure  6:  Parameter  estimates  as  a  function  of  the  number  of  frequencies.  Unmarked 
solid  lines  indicate  the  baseline  model  values.  Only  the  high  energy  part  of  the  signal 
is  utibzed  (from  the  250th  to  the  450th  sample). 
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Figure  7:  Parameter  estimates  as  a  function  of  the  number  of  frequencies.  Unmarked 
solid  lines  indicate  the  baseline  model  values.  The  time  series  is  diveded  into  five 
segments  and  an  averaging  method  is  employed. 
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Multi-frequency  matched  field 
processing  for  source  localization 


Based  on  the  same  data  set,  this  section  focuses  exclusively  on  the  localization 
performance  of  the  MFP  technique  for  three  reasons:  a)  in  this  experiment  the  source 
coordinates  are  known  and  thus  the  processor’s  performance  can  be  easily  verified,  b) 
the  prediction  of  the  source  location  can  be  used  as  an  indication  of  the  accuracy  with 
which  the  propagation  channel  is  simulated,  and  c)  to  demonstrate  the  significance 
of  the  accurate  knowledge  of  the  environment  in  localization  experiments. 

The  conclusions  derived  previously  are  employed  to  enhance  the  performance  of  the 
Bartlett-based  processor.  The  goal  is  to  find  the  most  probable  location  for  the  un¬ 
derwater  explosion  based  on  comparisons  between  the  signal  received  on  the  vertical 
array  and  synthetic  data.  Another  way  of  demonstrating  the  outcome  of  this  com¬ 
parison  is  by  plotting  ambiguity  surfaces  which  indicate  the  detection  performance 
for  all  possible  combinations  of  source  range  and  depth.  For  the  particular  problem, 
the  region  of  interest  is  from  11  to  13  km  in  range  and  from  10  to  20  m  in  depth. 
The  estimated  source  coordinates  are  at  a  range  of  12.1  km  and  depth  18  m. 

The  first  ambiguity  surface  shown  in  Fig  8a  represents  the  localization  performance 
for  the  single-frequency  case.  The  synthetic  data  are  modelled  according  to  the 
baseline  model.  The  environmental  parameters  are  assumed  to  be  fixed,  thus  the  in¬ 
version  is  confined  to  the  source  coordinates.  The  range  is  estimated  at  12.1  to  12.3 
km  while  the  depth  is  underestimated  at  11  to  16  m.  This  is  not  a  satisfactory  result 
since  the  area  of  uncertainty  remains  large.  To  improve  this  performance,  additional 
frequencies  are  incorporated  in  the  processor.  The  ten  best  frequencies  are  now  se¬ 
lected  because,  as  shown  in  section  5,  this  critical  number  yielded  stable  estimates 
for  the  particular  data  set.  Again,  the  synthetic  environment  remains  fixed  and  the 
comparison  between  real  and  synthetic  data  is  based  only  on  the  source  coordinates. 
Fig.  8b  shows  the  ambiguity  surface  for  this  case.  Although  the  estimated  range  is 
between  11.9  and  12.1  km,  there  is  no  depth  resolution.  The  fact  that  the  utiliza¬ 
tion  of  extra  frequencies  does  not  improve  the  processor’s  performance  implies  an 
inaccurately  estimated  synthetic  field  which  causes  the  additional  information  due 
to  combined  frequencies  to  have  a  misleading  effect  on  the  processor.  This  is  an 
indication  that  the  parameters  of  the  baseline  model  are  not  sufficiently  accurate. 
So,  the  next  step  is  to  generate  a  more  precise  representation  of  the  environment 
and  then  apply  the  same  multi-frequency  processor. 
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For  this  reason,  the  SAGA  algorithm  is  employed  as  a  pre-processor  that  estimates 
only  the  environmental  parameters.  The  set  which  provides  the  best  environmental 
lit  takes  the  place  of  the  baseline  model  and  the  ambiguity  surface  is  recalculated. 
Now  the  outcome  is  considerably  improved  compared  to  both  previous  cases.  As 
shown  in  Fig.  9a  the  most  probable  area  for  the  source  is  around  12.25  km  in  range 
and  between  16  and  19  m  in  depth. 

To  further  improve  the  processor’s  performance,  we  modified  the  sequential  opti¬ 
mization  procedure  to  a  global  and  simultaneous  one.  In  other  words,  instead  of  first 
calculating  the  best  possible  environment  and  then  searching  for  the  most  probable 
source  location,  now  the  inversion  is  jointly  performed  for  aU  environmental  and 
geometric  parameters  in  a  multi-frequency  mode.  The  outcome  is  a  five- dimensional 
a  posteriori  probability  function.  Figure  9b  shows  the  distributions  for  the  source 
coordinates.  The  only  peak  is  around  12.25  km  in  range  and  18  m  in  depth.  As  an¬ 
ticipated,  the  source  localization  performance  is  improved  by  using  multi-frequency 
global  inversion  algorithms. 
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Figure  8:  Ambiguity  surfaces  for  localization  performance  based  on  the  baseline 
model  (non-optimized  environment):  a)  single-frequency  case,  b)  multi-frequency 
case.  The  source  is  located  approximately  at  range  12.1  km  and  depth  18  m. 
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Figure  9:  Multi-frequency  localization  performance:  a)  ambiguity  surface  based  on 
the  optimized  environment,  b)  source  coordinate  distribution  surface  after  global 
inversion.  The  source  is  located  approximately  at  range  12.1  km  and  depth  18  m. 
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An  alternative  form  of  the  objective  function 


As  it  is  derived  in  section  3,  the  objective  function  is  the  summation  of  a  number 
of  frequency-dependent  terms.  Every  term  expresses  the  difference  between  the 
maximum  power  (perfect  match)  and  the  Bartlett  power  (actual  match).  Therefore, 
the  value  of  the  function  Eq.  (14)  varies  from  zero  to  the  sum  of  maximum  power. 
This  sum  is  an  unnormahzed  quantity  and  the  highest  value  it  can  assume  depends 
not  only  on  the  number  of  frequencies  involved  but  also  on  the  energy  content  of 
each  frequency.  As  the  energy  content  of  frequencies  varies  from  case  to  case,  direct 
comparison  between  different  scenarios  is  difficult  because  every  problem  has  its  own 
bounds. 

To  overcome  this  inconsistency,  a  normalized  version  of  the  objective  function  was 
previously  employed  [14].  As  seen  in  Eq.  (17), 


L 

/=! 

each  term  is  normalized  with  the  trace  of  the  covariance  matrix,  a  quantity  that 
represents  the  maximum  power  output  for  each  frequency.  Now  the  bounds  of  the 
objective  function  change  from  zero  (perfect  match)  to  the  number  of  frequencies 
used  (no  match  at  all).  Therefore,  different  scenarios  which  incorporate  the  same 
number  of  frequencies  can  be  compared  quantitatively  as  they  have  the  same  per¬ 
formance  limits. 

This  normaUzation  scheme  has  a  major  disadvantage:  it  eliminates  the  absolute 
difference  between  individual  frequency  components.  Since  every  Bartlett  term  is 
normalized  with  its  own  maximum  power,  frequency  components  which  encompass 
small  amounts  of  energy  are  represented  as  equal  to  those  belonging  to  the  high 
energy  part  of  the  spectrum.  Consequently,  the  parameter  estimates  can  be  very 
unstable,  as  demonstrated  in  Fig.  10.  Here,  the  source  coordinates  are  estimated 
using  both  the  normalized  Eq.  (14)  and  the  unnormalized  Eq.  (17)  version  of  the 
objective  function.  The  normalized-version  results  have  unpredictable  behaviour, 
especially,  in  the  case  of  source  depth  where  the  estimate  does  not  converge  at  aU. 
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Figure  10:  Parameter  estimates  for  source  range  (a)  and  for  source  depth  (b).  Solid 
lines  indicate  the  baseline  model  values.  Lines  designated  with  crosses  represent 
the  estimates  of  the  normalized  objective  function.  Lines  designated  with  circles 
indicate  the  estimates  of  the  unnormalized  objective  function. 


The  normalized  version  is  suitable  only  when  all  the  frequency  components  convey 
similar  amounts  of  energy,  and  thus  normalization  is  just  a  scaling  process  which 
does  not  affect  the  performance  of  the  processor.  When  this  condition  is  not  fulfilled, 
or  there  is  no  sufficient  a  priori  knowledge  about  the  frequency  components  to  be 
processed,  the  use  of  the  unnormalized  objective  function  is  recommended. 
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Conclusions 


In  processing  shot  data  in  shallow  water,  it  was  found  that  single  frequency  pa¬ 
rameter  estimates  are  unstable  and  unreliable.  The  problem  is  overcome  by  using 
multi-frequency  processing  schemes.  The  frequency  selection  is  based  on  high  en¬ 
ergy  areas  of  the  spectrum.  On  the  basis  of  a  limited  data  set,  it  is  shown  that  as 
the  number  of  frequencies  is  increased,  the  estimates  remain  unchanged  beyond  a 
critical  point.  This  observation  can  be  used  to  reduce  computation  time  because 
it  places  an  upper  limit  on  the  number  of  frequencies  to  be  incorporated  into  the 
processor.  The  “saturation”  point  of  ten  frequencies  found  here  is  expected  to  vary 
for  different  data  sets. 

It  is  demonstrated  that  the  localization  performance  of  the  MFP  techniques  is  en¬ 
hanced  by  using  multi-frequency  methods.  Unless  environmental  modelling  is  ac¬ 
curate,  the  utilization  of  many  frequencies  does  not  necessarily  imply  improved 
source  resolution.  The  combination  of  optimized  environmental  parameters  and 
multi-frequency  processing  techniques  enhances  the  performance  of  Bartlett-based 
processors. 

It  is  found  that  when  the  inversion  for  the  source  coordinates  is  incorporated  in  a 
global  search  scheme  which  also  allows  geoacoustic  parameters  to  be  jointly  opti¬ 
mized,  the  localization  performance  improves  considerably.  As  this  type  of  search  is 
extremely  time  consuming,  the  application  of  an  efficient  global  inversion  method, 
such  as  genetic  algorithms,  is  important. 
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